!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! The program that coded for solve the supersonic cone flow field 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Creator : B. G.
! Date    : 2015-12-05
! Revised : 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Input  : in.txt
! Output : out.txt
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Contain :
!   SHOCK
!   RK44
!   INTERSECTION
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Note:
!   I. all angle in rad, except for input and output
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

program TM
implicit none
! definition
integer::N,i,j,k
real*8::M1,p1,rho1,T1,V1,deltaC,gamma,pi,R,alpha1,M1A,M2A
real*8,allocatable::eps(:),psi(:),delPsi(:),uBA(:),vBA(:),deltaCC(:),u(:),v(:)
real*8,allocatable::theta(:),MCA(:),M(:),p(:),rho(:),T(:),Vi(:),x(:),y(:)
real*8::p2pt,r2rt,ppt,rrt

! input parameters
open(unit=101,file='in.txt',status='old')
read(101,*) N                     ! step number
read(101,*) M1                    ! freestream Mach number
read(101,*) p1                    ! freestream pressure
read(101,*) rho1                  ! freestream density
read(101,*) T1                    ! freestream temperature
read(101,*) deltaC                ! cone semiangle
allocate(x(N+1))                  ! streamline points (x)
allocate(y(N+1))                  ! streamline points (y)
read(101,*) y(1)                  ! starting point of streamline tracing (y)
close(101)                        ! close input file

! allocate array
allocate(eps(N))                ! shock wave angles 
allocate(psi(N+5))                ! angle steps between shock and cone surface
allocate(delPsi(N))               ! angle step values
allocate(uBA(2*N+5))              ! radical non-d critical velocity
allocate(vBA(2*N+5))              ! circumferential non-d critical velocity
allocate(deltaCC(N))              ! calculated cone semiangles
allocate(theta(N+1))              ! angle of velocity
allocate(u(N+1))                  ! velocity along x-direction (freestream)
allocate(v(N+1))                  ! velocity along y-direction
allocate(MCA(N+1))                ! critical Mach number
allocate(M(N+1))                  ! Mach number
allocate(p(N+1))                  ! pressure
allocate(rho(N+1))                ! density
allocate(T(N+1))                  ! temperature
allocate(Vi(N+1))                  ! velocity

! initialization
gamma=1.4                         ! specific heat ratio
pi=3.141592653
R=287.04                          ! gas constant, J/kg-K
deltaC=deltaC*pi/180              ! deltaC in rad
alpha1=asin(1/M1)                 ! freestream Mach angle
M1A=sqrt(((gamma+1)*M1*M1)/(2+(gamma-1)*M1*M1))
                                  ! freestream critical Mach number
V1=M1*sqrt(gamma*R*T1)            ! freestream velocity
                                  
do i=1,N                          ! i-loop start
   if (i==1) then                    ! i if start               
      eps(i)=deltaC+0.5*alpha1       ! shock wave angle assumption for first loop
   else if (i==2) then
      eps(i)=eps(i-1)+deltaC-deltaCC(i-1)
   else
      eps(i)=eps(i-1)+(eps(i-1)-eps(i-2))*(deltaC-deltaCC(i-1))/&
         &(deltaCC(i-1)-deltaCC(i-2))
                                     ! secant method for eps
   end if                            ! i if end 
   
   delPsi(1)=-(eps(i)-deltaC)/N      ! first psi step (negtive value)
   
   ! determine the flow properties immediately after the shock wave
   call SHOCK(M1,p1,rho1,T1,V1,M1A,eps(i),uBA(1),vBA(1),M2A,theta(1),p(1),rho(1),&
      &T(1),Vi(1))
   
   ! integrate to cone surface
   do j=1,N+5                        ! j-loop start
      psi(j)=eps(i)+(j-1)*delPsi(1)
      call RK44(uBA(j),vBA(j),delPsi(1),psi(j),uBA(j+1),vBA(j+1))
      ! vBA(j+1) >= 0 , exit the j-loop
      if ((abs(vBA(j+1)) <= 1e-6) .or. (vBA(j+1) > 1e-6)) exit
   end do                            ! j-loop end
   ! vBA(j+1) = 0  
   ! if (abs(vBA(j+1)) < 1e-6) exit
   deltaCC(i)=psi(j)+delPsi(1)          ! calculated cone semiangle
   psi(j+1)=psi(j)+delPsi(1)            ! calculated last psi
   ! vBA(j+1) > 0 
   if (abs(vBA(j+1)) > 1e-6) then         ! vBA(j+1) > 0, if start
      ! second, interpolation for delta_psi
      delPsi(2)=delPsi(1)*(0.0-vBA(j))/(vBA(j+1)-vBA(j))
      call RK44(uBA(j),vBA(j),delPsi(2),psi(j),uBA(j+2),vBA(j+2))
      deltaCC(i)=psi(j)+delPsi(2)
      psi(j+1)=psi(j)+delPsi(2)
      ! successive, secant method for delta_psi
      if (abs(vBA(j+2)) > 1e-6) then    ! abs(vBA(j+2)) > 1e-6, if start
         do k=j+3,j+3+N                    ! k-loop start
            delPsi(2+k-j-2)=delPsi(2+k-j-2-1)+(delPsi(2+k-j-2-1)-delPsi(2+k-j-2-2))*&
               &(0.0-vBA(k-1))/(vBA(k)-vBA(k-2))
            call RK44(uBA(j),vBA(j),delPsi(2+k-j-2),psi(j),uBA(k),vBA(k))
            deltaCC(i)=psi(j)+delPsi(2+k-j-2)
            psi(j+1)=psi(j)+delPsi(2+k-j-2)
            if (abs(vBA(k)) < 1e-6) exit
         end do                           ! k-loop end
      end if                           ! abs(vBA(j+2)) > 1e-6, if end
   end if                           ! vBA(j+1) > 0, if end
   
   ! determine if deltaCC(i) equals to deltaC
   if (abs(deltaCC(i)-deltaC) <= 1e-6) exit
   
end do                           ! i-loop end

! determine flow properties within the flowfield calculated
u(1)=Vi(1)*cos(theta(1))
v(1)=Vi(1)*sin(theta(1))
M(1)=sqrt((2*M2A*M2A/(gamma+1))/(1-(gamma-1)*M2A*M2A/(gamma+1)))
MCA(1)=M2A
x(1)=y(1)/tan(psi(1))
p2pt=(1-(gamma-1)*MCA(1)*MCA(1)/(gamma+1))**(gamma/(gamma-1))
r2rt=(1-(gamma-1)*MCA(1)*MCA(1)/(gamma+1))**(1/(gamma-1))
do j=2,N+1                      ! j-loop start
   MCA(j)=sqrt(uBA(j)*uBA(j)+vBA(j)*vBA(j))
   M(j)=sqrt((2*MCA(j)*MCA(j)/(gamma+1))/(1-(gamma-1)*MCA(j)*MCA(j)/(gamma+1)))
   ppt=(1-(gamma-1)*MCA(j)*MCA(j)/(gamma+1))**(gamma/(gamma-1))
   rrt=(1-(gamma-1)*MCA(j)*MCA(j)/(gamma+1))**(1/(gamma-1))
   p(j)=p(1)*ppt/p2pt
   rho(j)=rho(1)*rrt/r2rt
   T(j)=T(1)*p(j)*rho(1)/(p(1)*rho(j))
   Vi(j)=M(j)*sqrt(gamma*R*T(j))
   theta(j)=psi(j)-atan(abs(vBA(j))/uBA(j))
   call INTERSECTION(x(j-1),y(j-1),theta(j-1),psi(j),x(j),y(j))
                               ! calculate streamline point
   u(j)=Vi(j)*cos(theta(j))
   v(j)=Vi(j)*sin(theta(j))
end do                         ! j-loop end

! output
open(unit=201,file='output.txt')
write(201,2001)
   2001 format(/,T10,'Taylor-Maccoll flow around a cone'//,T4,'M1',16x,'p1 / Pa',&
      &11x,'rho1 / kg/m^3',5x,'T1 / K',12x,'deltaC / deg')
write(201,*) M1,p1,rho1,T1,deltaC*180/pi
write(201,*)
write(201,2002)
   2002 format(T4,'i',2x,'psi',15x,'M',17x,'p',17x,'rho',15x,'T',17x,'V',17x,&
      &'u',17x,'v',17x,'x',17x,'y')
do j=1,N+1                     ! j-loop start
   write(201,*) j,psi(j)*180/pi,M(j),p(j),rho(j),T(j),Vi(j),u(j),v(j),x(j),y(j)
end do                         ! j-loop end
close(201)

end program TM
